﻿using System;
using LatoolNet;

namespace NumericalSolution {
  class Poisson2DFVIterative : Poisson2DFV {

    private double[,] m_u;

    public override void Prepare() {
      m_u = new double[p_numberOfControlVolumes, p_numberOfControlVolumes];
      for (int i = 0; i < p_numberOfControlVolumes; i++) {
        for (int j = 0; j < p_numberOfControlVolumes; j++) {
          m_u[i, j] = 0.0;
        }
      }
    }

    public override ROResult2D Solve() {

      for (int i = 0; i < p_numberOfControlVolumes; i++) {

        double[] cprime = new double[p_numberOfControlVolumes];
        double[] dprime = new double[p_numberOfControlVolumes];

        for (int j = 0; j < p_numberOfControlVolumes; j++) {

          double coef_w = Double.NaN;
          double coef_e = Double.NaN;
          double coef_n = Double.NaN;
          double coef_s = Double.NaN;
          double coef_p = 0.0;

          if (i == 0) {
            coef_p += 2 * p_deltay / p_deltax;
            coef_w = 0.0;
          } else {
            coef_p += p_deltay / p_deltax;
            coef_w = p_deltay / p_deltax * m_u[i - 1, j];
          }

          if (i == p_numberOfControlVolumes - 1) {
            coef_p += 2 * p_deltay / p_deltax;
            coef_e = 0.0;
          } else {
            coef_p += p_deltay / p_deltax;
            coef_e = p_deltay / p_deltax * m_u[i + 1, j];
          }

          if (j == 0) {
            coef_p += 2 * p_deltax / p_deltay;
            coef_s = 0.0;            
          } else {
            coef_p += p_deltax / p_deltay;
            coef_s = p_deltax / p_deltay;
          }

          if (j == p_numberOfControlVolumes - 1) {
            coef_p += 2 * p_deltax / p_deltay;
            coef_n = 0.0;
          } else {
            coef_p += p_deltax / p_deltay;
            coef_n = p_deltax / p_deltay;
          }

          double a = -coef_s;
          double b = coef_p;
          double c = -coef_n;
          double d = coef_w +
                     coef_e +
                     f(posx(i), posy(j)) * p_deltax * p_deltay;

          if (j == 0) {
            cprime[j] = c / b;
          } else {
            cprime[j] = c / (b - a * cprime[j - 1]);
          }

          if (j == 0) {
            dprime[j] = d / b;
          } else {
            dprime[j] = (d - a * dprime[j - 1]) / (b - a * cprime[j - 1]);
          } 
        }

        m_u[i, p_numberOfControlVolumes - 1] 
          = dprime[p_numberOfControlVolumes - 1];
        for (int j = p_numberOfControlVolumes - 2; j >= 0; j--) {
          m_u[i, j] =  dprime[j] - cprime[j] * m_u[i, j + 1];
        }
      }
      return new ROResult2D(m_u);
    }
  }
}
